home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / examples / misc / tornado.pro < prev    next >
Encoding:
Text File  |  1997-07-08  |  12.6 KB  |  354 lines
  1. PRO TORNADO
  2.  
  3.    winx = 640.0
  4.    winy = 512.0
  5.  
  6.    clouds_pts = 200
  7.    ground_pts = 200
  8.    tornado_pts = 25
  9.    arc_pts = 7
  10.  
  11.    clouds_center_x = 0.5
  12.    clouds_center_y = 0.5
  13.  
  14.    trans = Fltarr(4, 4)
  15.    ident4 = trans
  16.    ident4([0, 5, 10, 15]) = 1.0
  17.  
  18.    max_torn_pt = tornado_pts - 1
  19.    max_torn_pt_m1 = max_torn_pt - 1
  20.    max_torn_pt_m2 = max_torn_pt - 2
  21.  
  22.    arc_ones = Replicate(1.0, arc_pts)
  23.    zeroes = Fltarr(tornado_pts)
  24.    ramp = 1.0 - (Findgen(tornado_pts) / Float(tornado_pts - 1))
  25.    clouds_ang = Fltarr(clouds_pts)
  26.    ground_ang = Fltarr(ground_pts)
  27.  
  28.    temp = Randomu(s, 1)
  29.  
  30.    pos_clouds_x = (Randomu(s, clouds_pts) * 2.0) - 0.5
  31.    pos_clouds_y = (Randomu(s, clouds_pts) * 2.0) - 0.5
  32.    pos_clouds_z = Replicate(1.0, clouds_pts)
  33.  
  34.    vel_clouds_x = Randomu(s, clouds_pts) /  100.0
  35.    vel_clouds_y = Randomu(s, clouds_pts) /  100.0
  36.    vel_clouds_z = Randomu(s, clouds_pts) / 1000.0
  37.  
  38.    pos_ground_x = (Randomu(s, ground_pts) * 2.0) - 0.5
  39.    pos_ground_y = (Randomu(s, ground_pts) * 2.0) - 0.5
  40.    pos_ground_z = Fltarr(ground_pts)
  41.  
  42.    vel_ground_x = Fltarr(ground_pts)
  43.    vel_ground_y = Fltarr(ground_pts)
  44.    vel_ground_z = Fltarr(ground_pts)
  45.  
  46.    pos_tornado_x = Replicate(clouds_center_x, tornado_pts)
  47.    pos_tornado_y = Replicate(clouds_center_y, tornado_pts)
  48.    pos_tornado_z = Replicate(1.0, tornado_pts)
  49.  
  50.    vel_tornado_x = Fltarr(tornado_pts)
  51.    vel_tornado_y = Fltarr(tornado_pts)
  52.    vel_tornado_z = Replicate((-0.01), tornado_pts) * ramp
  53.  
  54.    twist_ang = Randomu(s, tornado_pts) * 180.0
  55.    ang_inc = 360.0 / Float(arc_pts)
  56.  
  57.    spin_pts = Fltarr(3, arc_pts, tornado_pts)
  58.  
  59.    plane_pts = Fltarr(3, arc_pts)
  60.    ang = 0.0
  61.    FOR i=0, (arc_pts - 1) DO BEGIN
  62.       plane_pts(0, i) = Cos(ang * !DTOR)
  63.       plane_pts(1, i) = Sin(ang * !DTOR)
  64.       ang = ang + ang_inc
  65.    ENDFOR
  66.  
  67.    spin_tornado = Fltarr(tornado_pts)
  68.    spin_change = Replicate(0.1, tornado_pts)
  69.    speed_tornado = Fltarr(tornado_pts)
  70.  
  71.    Window, /Free, /Pixmap, Xsize=winx, Ysize=winy, Retain=1
  72.    pix_win = !D.Window
  73.  
  74.    Window, /Free, Title="TORNADO !", Xsize=winx, Ysize=winy, Retain=1
  75.    see_win = !D.Window
  76.    c_back = 0
  77.    c_ground = 1
  78.    c_clouds = 2
  79.    c_torn = 3
  80.    c_spot = 4
  81.    Tvlct, 015, 015, 031, c_back
  82.    Tvlct, 191, 063, 031, c_ground
  83.    Tvlct, 063, 063, 255, c_clouds
  84.    Tvlct, 255, 255, 255, c_torn
  85.    Tvlct, 255, 000, 000, c_spot
  86.  
  87.    Create_View, Winx=winx, Winy=winy, Az=(30), Ay=(0), Ax=(-90), $
  88.                 Zfac=1.0, Zoom=0.6, Persp=2.0
  89.  
  90.    time_inc = 0.025
  91.    float_vel = 0.1
  92.    grav = 1.0
  93.    grav_time = grav * time_inc
  94.    pi2 = !PI * 2.0
  95.  
  96.    count = 0.0
  97.    WHILE (1) DO BEGIN
  98.  
  99.       count = count + 1.0
  100.  
  101.       vel_clouds_z = Temporary(vel_clouds_z) + (float_vel / 20.0)
  102.  
  103.       z_pos_ramp = 1.0 - pos_tornado_z
  104.  
  105.       vel_tornado_x = Temporary(vel_tornado_x) + $
  106.          (ramp * z_pos_ramp^2 * (Randomu(s, tornado_pts) - 0.5) / 4.0)
  107.  
  108.       vel_tornado_y = Temporary(vel_tornado_y) + $
  109.          (ramp * z_pos_ramp^2 * (Randomu(s, tornado_pts) - 0.5) / 4.0)
  110.  
  111.       vel_tornado_z = Temporary(vel_tornado_z) - $
  112.          (ramp * float_vel * (Randomu(s, tornado_pts) - 0.5))
  113.  
  114.       IF (pos_tornado_z(0) GE 0.25) THEN $
  115.          vel_tornado_z(0) = vel_tornado_z(0) < vel_tornado_z(1) < (0.0)
  116.  
  117.       vel_tornado_x(max_torn_pt) = 0.0
  118.       vel_tornado_y(max_torn_pt) = 0.0
  119.       vel_tornado_z(max_torn_pt) = 0.0
  120.  
  121.       vel_tornado_x = Smooth(Temporary(vel_tornado_x),  3, /Edge_Truncate)
  122.       vel_tornado_y = Smooth(Temporary(vel_tornado_y),  3, /Edge_Truncate)
  123.       vel_tornado_z = Smooth(Temporary(vel_tornado_z),  3, /Edge_Truncate)
  124.  
  125.       pos_tornado_x = Temporary(pos_tornado_x) + (vel_tornado_x * time_inc)
  126.       pos_tornado_y = Temporary(pos_tornado_y) + (vel_tornado_y * time_inc)
  127.       pos_tornado_z = Temporary(pos_tornado_z) + (vel_tornado_z * time_inc)
  128.  
  129. ;     IF (pos_tornado_z(0) GT 0.0) THEN BEGIN
  130. ;        index = Where((pos_tornado_z(1:*)-pos_tornado_z(0:max_torn_pt_m1)) LE 0.0)
  131. ;        IF (index(0) GT 0L) THEN $
  132. ;           pos_tornado_z(index) = pos_tornado_z(index) - (0.005)
  133. ;     ENDIF
  134.  
  135.       pos_tornado_x(0) = pos_tornado_x(1)
  136.       pos_tornado_x(max_torn_pt) = pos_tornado_x(max_torn_pt_m1)
  137.       pos_tornado_x = Smooth(Temporary(pos_tornado_x),  3, /Edge_Truncate)
  138.  
  139.       pos_tornado_y(0) = pos_tornado_y(1)
  140.       pos_tornado_y(max_torn_pt) = pos_tornado_y(max_torn_pt_m1)
  141.       pos_tornado_y = Smooth(Temporary(pos_tornado_y),  3, /Edge_Truncate)
  142.  
  143.       pos_tornado_z(0) = pos_tornado_z(0) < (pos_tornado_z(1) - 0.01)
  144.       pos_tornado_z(max_torn_pt) = 1.0
  145.       pos_tornado_z = Smooth(Temporary(pos_tornado_z),  3)
  146.  
  147.       index = Where((pos_tornado_x LE 0.0) OR (pos_tornado_x GE 1.0))
  148.       IF (index(0) GE 0L) THEN vel_tornado_x(index) = 0.0
  149.       index = Where((pos_tornado_y LE 0.0) OR (pos_tornado_y GE 1.0))
  150.       IF (index(0) GE 0L) THEN vel_tornado_y(index) = 0.0
  151.       index = Where((pos_tornado_z LE 0.0) OR (pos_tornado_z GE 1.0))
  152.       IF (index(0) GE 0L) THEN vel_tornado_z(index) = 0.0
  153.       pos_tornado_x = (Temporary(pos_tornado_x) > 0.0) < 1.0
  154.       pos_tornado_y = (Temporary(pos_tornado_y) > 0.0) < 1.0
  155.       pos_tornado_z = (Temporary(pos_tornado_z) > 0.0) < 1.0
  156.  
  157. ;     speed_tornado =  ((1.5 + ramp) * time_inc * 20.0 * (Sqrt(spin_tornado) + 0.001))
  158.       speed_tornado =  ((0.5 + ramp) * time_inc * 50.0 * (Sqrt(spin_tornado) + 0.001))
  159.  
  160.       spin_change = Temporary(spin_change) + $
  161.          Smooth(((Randomu(s, tornado_pts) - 0.5) / 10.0), 3, /Edge_Truncate)
  162.       spin_change = Smooth(Temporary(spin_change), 3, /Edge_Truncate)
  163.  
  164.       spin_change = (Temporary(spin_change) > (-0.125)) < (0.125)
  165.       spin_tornado = ((Temporary(spin_tornado) + (spin_change * time_inc)) > 0.0) < 0.15
  166.       spin_tornado(max_torn_pt) = $
  167.          spin_tornado(max_torn_pt_m1) + ((0.001 * count) < 0.15)
  168.  
  169.       spin_tornado = Smooth(Temporary(spin_tornado),  3, /Edge_Truncate)
  170.       spin_tornado = Smooth(Temporary(spin_tornado), (tornado_pts/4), /Edge_Truncate)
  171.       spin_tornado(0) = spin_tornado(1) * (1.0 - pos_tornado_z(0))^4
  172.  
  173.       twist_ang = Temporary(twist_ang) + (10.0 * time_inc * speed_tornado)
  174.       index = Where(twist_ang GE 360.0)
  175.       IF (index(0) GE 0L) THEN twist_ang(index) = twist_ang(index) - 360.0
  176.  
  177.       ground_dx = pos_ground_x - pos_tornado_x(0)
  178.       ground_dy = pos_ground_y - pos_tornado_y(0)
  179.       ground_dz = pos_ground_z - pos_tornado_z(0)
  180.       ground_dist = Sqrt(ground_dx^2 + ground_dy^2)
  181.       index = Where(ground_dist GT 0.0)
  182.       ground_ang(*) = 0.0
  183.       ground_ang(index) = Atan(ground_dy(index), ground_dx(index))
  184.  
  185.       ground_xyz_dist = Sqrt(ground_dx^2 + ground_dy^2 + ground_dz^2)
  186.  
  187.       vect_x = pos_tornado_x(2) - pos_tornado_x(0)
  188.       vect_y = pos_tornado_y(2) - pos_tornado_y(0)
  189.       pt_dz = Sqrt(vect_x^2 + vect_y^2)
  190.       IF (pt_dz GT 0.0) THEN pt_ang = Atan(vect_y, vect_x) ELSE pt_ang = 0.0
  191.       ang_diff = pt_ang - ground_ang
  192.       tilt_fac = (pt_dz * (-Cos(ang_diff))) + (Randomu(s, ground_pts) / 2.0)
  193.  
  194.       vel_ground_z = Temporary(vel_ground_z) - (grav_time)
  195.  
  196.       pos_ground_x = ((Temporary(pos_ground_x) + (vel_ground_x * time_inc)) > (-1.0)) < 2.0
  197.       pos_ground_y = ((Temporary(pos_ground_y) + (vel_ground_y * time_inc)) > (-1.0)) < 2.0
  198.       pos_ground_z = (Temporary(pos_ground_z) + (vel_ground_z * time_inc)) > 0.0
  199.       index = Where(pos_ground_z LE 0.0)
  200.       IF (index(0) GE 0L) THEN BEGIN
  201.          vel_ground_x(index) = 0.0
  202.          vel_ground_y(index) = 0.0
  203.          vel_ground_z(index) = time_inc * spin_tornado(1) * 10.0 * tilt_fac(index) / $
  204.             (ground_xyz_dist(index) + 0.25)^6
  205.       ENDIF
  206.  
  207.       index = Where(pos_ground_z GT 0.0)
  208.       IF (index(0) GE 0L) THEN BEGIN
  209.          vel_ground_x(index) = vel_ground_x(index) - $
  210.             (Sin(ground_ang(index)) * 4.0 * spin_tornado(1) / $
  211.             (1.5*(ground_xyz_dist(index) + 0.5))^6)
  212.          vel_ground_y(index) = vel_ground_y(index) + $
  213.             (Cos(ground_ang(index)) * 4.0 * spin_tornado(1) / $
  214.             (1.5*(ground_xyz_dist(index) + 0.5))^6)
  215.       ENDIF
  216.  
  217.       vel_ground_x = Temporary(vel_ground_x) * 0.85
  218.       vel_ground_y = Temporary(vel_ground_y) * 0.85
  219.       vel_ground_z = Temporary(vel_ground_z) * 0.95
  220.  
  221.       clouds_dx = pos_clouds_x - pos_tornado_x(max_torn_pt)
  222.       clouds_dy = pos_clouds_y - pos_tornado_y(max_torn_pt)
  223.       clouds_dist = Sqrt(clouds_dx^2 + clouds_dy^2)
  224.       index = Where(clouds_dist GT 0.0)
  225.       clouds_ang(*) = 0.0
  226.       clouds_ang(index) = Atan(clouds_dy(index), clouds_dx(index))
  227.       vel_clouds_x = Temporary(vel_clouds_x) - $
  228.          (Sin(clouds_ang) * 2.0 * spin_tornado(max_torn_pt) / (1.5*(clouds_dist + 0.5))^4)
  229.       vel_clouds_y = Temporary(vel_clouds_y) + $
  230.          (Cos(clouds_ang) * 2.0 * spin_tornado(max_torn_pt) / (1.5*(clouds_dist + 0.5))^4)
  231.       vel_clouds_x = Temporary(vel_clouds_x) + (0.075 * (Randomu(s, clouds_pts) - 0.5))
  232.       vel_clouds_y = Temporary(vel_clouds_y) + (0.075 * (Randomu(s, clouds_pts) - 0.5))
  233.       vel_clouds_x = Temporary(vel_clouds_x) * 0.75
  234.       vel_clouds_y = Temporary(vel_clouds_y) * 0.75
  235.  
  236.       vel_clouds_x = Temporary(vel_clouds_x) - $
  237.          (0.03 * Cos(clouds_ang) / (1.0 + (10.0 * clouds_dist)))
  238.       vel_clouds_y = Temporary(vel_clouds_y) - $
  239.          (0.03 * Sin(clouds_ang) / (1.0 + (10.0 * clouds_dist)))
  240.  
  241.       vect_x = pos_tornado_x(max_torn_pt) - pos_tornado_x(max_torn_pt_m2)
  242.       vect_y = pos_tornado_y(max_torn_pt) - pos_tornado_y(max_torn_pt_m2)
  243.       pt_dz = Sqrt(vect_x^2 + vect_y^2)
  244.       IF (pt_dz GT 0.0) THEN pt_ang = Atan(vect_y, vect_x) ELSE pt_ang = 0.0
  245.       ang_diff = pt_ang - clouds_ang
  246.       tilt_fac = pt_dz * (-Cos(ang_diff)) * 5.0
  247.  
  248.       clouds_dz = pos_clouds_z - pos_tornado_z(max_torn_pt)
  249.       clouds_dist = (clouds_dx^2 + clouds_dy^2 + clouds_dz^2)
  250.  
  251.       vel_clouds_z = Temporary(vel_clouds_z) - $
  252.          (spin_tornado(max_torn_pt) * tilt_fac / (clouds_dist > 0.01))
  253.  
  254.       vel_clouds_z = Temporary(vel_clouds_z) + $
  255.          ((Randomu(s, clouds_pts) - 0.5) * float_vel / 40.0)
  256.  
  257.       vel_clouds_z = (Temporary(vel_clouds_z) > (-float_vel)) < float_vel
  258.  
  259.       pos_clouds_x = Temporary(pos_clouds_x) + (vel_clouds_x * time_inc)
  260.       pos_clouds_y = Temporary(pos_clouds_y) + (vel_clouds_y * time_inc)
  261.       pos_clouds_z = ((Temporary(pos_clouds_z) + (vel_clouds_z * time_inc)) > 0.0) < 1.0
  262.  
  263.       index = Where(Sqrt((pos_clouds_x - 0.5)^2 + (pos_clouds_y - 0.5)^2) GT 1.0)
  264.       IF (index(0) GE 0L) THEN BEGIN
  265.         vel_clouds_x(index) = 0.0
  266.         vel_clouds_y(index) = 0.0
  267.       ENDIF
  268.  
  269.  
  270.       Wset, pix_win
  271.       Erase, c_back
  272.  
  273.       FOR i=0, max_torn_pt DO BEGIN
  274.  
  275.          prev_pt = (i - 1) > 0
  276.          next_pt = (i + 1) < max_torn_pt
  277.  
  278.          vect_x = pos_tornado_x(next_pt) - pos_tornado_x(prev_pt)
  279.          vect_y = pos_tornado_y(next_pt) - pos_tornado_y(prev_pt)
  280.          vect_z = (pos_tornado_z(next_pt) - pos_tornado_z(prev_pt)) > 0.001 > $
  281.                                            (0.01 * Float(i) / Float(max_torn_pt))
  282.          xy_len = Sqrt(vect_x^2 + vect_y^2)
  283.  
  284.          IF (xy_len GT 0.0) THEN BEGIN
  285.             ang_z = Atan(vect_y, vect_x)
  286.             ang_y = (!PI / 2.0) - (Atan(vect_z, xy_len))
  287.          ENDIF ELSE BEGIN
  288.             ang_z = 0.0
  289.             ang_y = 0.0
  290.          ENDELSE
  291.  
  292.          trans = ident4
  293.  
  294.          sf = (spin_tornado(i)) * 1.25
  295.          m4x4 = ident4
  296.          m4x4([0, 5, 10]) = [sf, sf, 1.0]
  297.          trans = Temporary(trans) # m4x4
  298.  
  299.          m4x4 = ident4
  300.          s_ang = Sin(twist_ang(i)-ang_z)
  301.          c_ang = Cos(twist_ang(i)-ang_z)
  302.          m4x4(0,0) = c_ang
  303.          m4x4(0,1) = s_ang
  304.          m4x4(1,0) = (-s_ang)
  305.          m4x4(1,1) = c_ang
  306.          trans = Temporary(trans) # m4x4
  307.  
  308.          m4x4 = ident4
  309.          s_ang = Sin(ang_y)
  310.          c_ang = Cos(ang_y)
  311.          m4x4(0,0) = c_ang
  312.          m4x4(0,2) = (-s_ang)
  313.          m4x4(2,0) = s_ang
  314.          m4x4(2,2) = c_ang
  315.          trans = Temporary(trans) # m4x4
  316.  
  317.          m4x4 = ident4
  318.          s_ang = Sin(ang_z)
  319.          c_ang = Cos(ang_z)
  320.          m4x4(0,0) = c_ang
  321.          m4x4(0,1) = s_ang
  322.          m4x4(1,0) = (-s_ang)
  323.          m4x4(1,1) = c_ang
  324.          trans = Temporary(trans) # m4x4
  325.  
  326.          m4x4 = ident4
  327.          m4x4([3,7,11]) = [pos_tornado_x(i), pos_tornado_y(i), pos_tornado_z(i)]
  328.          trans = Temporary(trans) # m4x4
  329.  
  330.          new_pts = Vert_T3d(plane_pts, Matrix=trans)
  331.          new_pts(2, *) = (new_pts(2, *) > 0.0) < 1.0
  332.          spin_pts(0, 0, i) = new_pts
  333.  
  334.       ENDFOR
  335.  
  336.       Plots, pos_ground_x, pos_ground_y, pos_ground_z, Psym=4, Symsize=1.5, $
  337.          /Normal, /T3d, Color=c_ground
  338.  
  339.       Plots, pos_clouds_x, pos_clouds_y, pos_clouds_z, Psym=1, Symsize=1.0, $
  340.          /Normal, /T3d, Color=c_clouds
  341.  
  342.       Plots, Reform(spin_pts, 3, arc_pts*tornado_pts), $
  343.          Psym=4, Symsize=0.5, /Normal, /T3d, Color=c_torn
  344.  
  345. ;     Plots, pos_tornado_x, pos_tornado_y, pos_tornado_z, $
  346. ;        /Normal, /T3d, Color=c_spot
  347.  
  348.       Wset, see_win
  349.       Device, Copy=[0, 0, winx, winy, 0, 0, pix_win]
  350.       Empty
  351.    ENDWHILE
  352.  
  353. END
  354.